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Appendix: Successive Fourier Analysis A.l 
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Nomenclature 





Hij 


j 

k 

M 

n 


N 


PD: 


p « 


coefficient of cosine Fourier series 
coefficient of sine Fourier series 

the average value of the constant terms in the harmonic 

oscillation responses 

coefficient of a 1 term in static flow 

2- D lift coefficient. 

3- D lift coefficient. 

variation of lift coefficient with respect to angle of attack 
variation of lift coefficient with respect to pitch rate 
a constant 

constants associated with the virtual mass effect 
constants in amplitude function to be determined 
imaginary part of a complex number 
index 

reduced frequency (=co J/vJ 
Mach number 

index for reduced frequency, also index for the coefficient 

in Pad6 approximants 

the number of frequency 

Pade approximants 

coefficients for Pad6 approximants 


t 


time 


t’ 


UQVLM 

V oo 

Greek: 

a 

a l 

a 0 

a m 

a 

e 

x 

\ 

0 


nondimensional time (stv^l) 

unsteady quasi-vortex lattice method program 

free stream velocity 

angle of attack (=ocq coskt’) 

defined as c^+a 

amplitude of angle of attack 

mean angle of attack 

time rate of change in angle of attack 

reference length 

dummy time integration variable 
running variable in time 
defined as 0=kt’ 



Chapter 1 
Introduction 


Due to the requirement of increased performance and maneuverability, 
the flight envelope of a modem fighter is frequently extended to the high 
angle-of-attack regime. Vehicles maneuvering in this regime are subjected 
to nonlinear aerodynamic loads. The nonlinearities are due mainly to 
three-dimensional separated flow and concentrated vortex flow that occur 
at large angles of attack. Accurate prediction of these nonlinear airloads 
is of great importance in the analysis of a vehicle’s flight motion and in the 
design of its flight control system. As Tobak and Schiff mentioned in ref. 
1, the main difficulty in determining the relationship between the 
instantaneous aerodynamic load on a maneuvering vehicle and the motion 
variables is that this relationship is determined not only by the 
instantaneous values of motion variables but also by all of the prior states 
of the motion up to the current state. Due to advanced computer 
techniques, one straightforward way is to solve the flow-field problem and 
the dynamic equation together. For example, a CFD method can be used 
to solve the Navier-Stokes equations governing the separated flow field. 
Then the calculated forces and moments are used in the dynamic equations 
governing the vehicle’s motion to calculate motion variables. The motion 
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variables will change the vehicle’s attitude, and thus the forces and 
moments. Results of repeatedly calculating these coupled equations would 
be the complete time histories of the aerodynamic response and of the 
vehicle’s motion. Although solving these coupled equations is the exact way 
to account for the time-history effects in predicting the aerodynamic 
response to arbitrary maneuvers, this is obviously a very costly approach. 
In particular, at high angles of attack, the aerodynamic loads depend 
nonlinearly on the motion variables. Under such conditions, even if the 
vehicles start from closely similar initial conditions, they can experience 
widely varying motion histories. Thus, a satisfactory evaluation of the 
performance envelope of the aircraft may require a large number of coupled 
computations, one for each change in initial conditions. Further, since the 
motion and the aerodynamic response are linked together in this approach, 
there can be no reutilization of the previously obtained aerodynamic 
reactions. 

To avoid the disadvantage of solving the coupled flow- field equations and 
aircraft’s motion equations, an alternate approach is to use a mathematical 
modeling to describe the steady and unsteady aerodynamics for the 
aircraft’s equations of motion. Ideally, with a mathematical model, an 
evaluation of the aerodynamic terms specified by the model would be 
required only once. The specified model can be reutilized to solve the 
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aircraft’s equations of motion over a range of motion variables and flight 
conditions. 

In the classical linear potential flow theory (refs. 2 and 3), researchers 
in the field of aeroelasticity used the Fourier transform to relate the 
aerodynamic response of step change in angle of attack of a wing to that of 
harmonic oscillatory motions. The transient aerodynamic reaction to a step 
change is termed the "indicial function" and has been calculated for several 
classes of isolated wings (refs. 2-5). By a suitable superposition (ref. 6) of 
these results, the aerodynamic forces and moments induced in any 
maneuvers can be studied (refs. 2 and 3). Tobak has applied the indicial 
function concept to analyze the motions of wings and wing-tail combinations 
(ref. 7). Later, based on a consideration of function, Tobak and his 
colleagues (refs. 1 and 8) have extended the concept of indicial function into 
the nonlinear aerodynamic regimes. The simplest nonlinear aerodynamic 
model proposed in ref.l has been applied by several authors (refs. 9-13) to 
perform the analysis. However, that simplest model is accurate only to the 
first order of frequency. It needs to be improved for a more general 
response. 

Aerodynamic forces and moments acting on a rapidly maneuvering 
aircraft are, in general, nonlinear functions of motion variables, their time 
rate of change, and the history of maneuvering. How these unsteady 
aerodynamic forces and moments may be represented becomes uncertain, 
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in particular at high angles of attack. If the response is measured by wind- 
tunnel dynamic testing, questions arise as to how the measured time- 
history data can be analyzed and expressed in a form suitable for flight 
dynamic simulation. For a certain type of nonlinearities produced in a test 
with small-amplitude oscillation, the analysis has been accomplished by 
separating the time-history data into in-phase and out-phase components 
(ref. 14). When large-amplitude forced oscillations are employed in the 
wind- tunnel testing at a large mean angle of attack, the aerodynamic 
phenomena may involve dynamic stall and/or strong vortex flow, with or 
without vortex breakdown. In this case, higher harmonic components in 
the aerodynamic response are expected to exist (ref. 15) and the 
phenomenon of aerodynamic lag may be important. Therefore, a more 
general modeling technique is needed. 

In this research, a numerical method will be developed to analyze the 
nonlinear and time-dependent aerodynamic response to establish the 
generalized indicial function in terms of motion variables and their time 
rates of change. 
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Chapter 2 

Mathematical Development 


2.1 Aerodynamic Modeling 

Tobak and SchifF (ref. 1), based on a consideration of functional, 
developed a fundamental formulation of aerodynamic response. For 
example, the lift response to pitching oscillation may be given in the form 
of a generalized indicial response as 


C L (fc)=C L (0)+f C C^[t,T;a(5) ,d(5) ,<?(£) ,g(5)l -^-dx 


— f C c L [t,t;a(5),4 (5), <7(5). $(5)1-^* (l) 

v. Jo a ax 


where a and q are the time rate of change in angle of attack and the pitch 
rate, respectively, and t is the time. V is the free stream velocity and fi is 
a reference length. The variable ^ is a running variable in time over the 
interval 0 to x. This means that the indicial response depends not only on 
the current values of motion variables, but also on the past history of these 
variables. For practical implementation, eq. (1) requires further 
simplification. By introducing the assumption of a slowly varying motion, 
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Tobak and Schiff neglected the dependence of the indicial response on 
a and q. A slightly simplified expression of eq. (1) can be written as 


C L (t)=C L (0) 


+f o C C L '[t-x:a (t) a (x) ,q( t) ] da J^ ^ dx 

+ -nrf C c L [t-x;a(x)&(x) ,q(x )] d ^L dx 
V m J o a ax 


( 2 ) 


Although the form of eq. (2) represents a great simplification over that 
of eq. (1), the equation still includes the full linear form as a special case. 

The main objective in the present investigation is to find a suitable form 
for the integrand of eq. (2). Then the time response C L (t) can be calculated 
through the integration of eq. (2) by substituting the suitable form of 

cx 

and C Lq . In wind-tunnel testing, q is the same as a . Since the method 
developed in this study will be used to analyze wind tunnel data, a will be 
used instead of q in the following investigation and the investigation will 
be focused on lift force. 

In the linear theory (refs. 2 and 3) , the aerodynamic response could be 
separated into a product of an amplitude function and a phase function in 
harmonic motion. The amplitude function depends on motion variables and 
their time rate of change. On the other hand, the phase function is a 
function of frequency and accounts for any phase lag between the response 
and the excitation. In a two-dimensional linear theory, the phase function 
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is given by Theodorsen’s circulation function (refs. 2 and 3). After response 
has been obtained at different frequencies with the same amplitude in 
harmonic oscillation, the phase function can be determined numerically. 
After use of reciprocal relations (ref. 16), the indicial function can be 
defined by numerical means. This approach has been used for numerical 
determination of indicial lift for plunging airfoil in ref. 5 and for plunging 
wings in ref. 17. 

The method for the linear theory will be generalized as follows. 
Instead of assuming that the aerodynamic response is a product of an 
amplitude function and a phase function as it is in the linear theory, it is 
taken to be a sum of the products of amplitude functions and phase 
functions in harmonic motion; i.e., 


C,= C 0+ £ (.amplitude function) j * ( phase function) j 

In the linear theory, j equals 1 in the equation. To determine what is 
the form of the amplitude functions and the phase functions, the 
aerodynamic response due to harmonic oscillation is assumed to be of the 
form 


C L = F 0 +F(a,a )a + G(a,& )<x 


(3) 


7 



and it is defined that 


a l = + a o cos(kt’) 

a = a 0 cos(kt’) 
ex = (-a 0 k) sin(kt’) 

where k is the reduced frequency, t’ is the nondimensionalized time, is 
the mean angle of attack and a Q is the amplitude of angle of attack. To 
find the constant Fq and functions F and G as functions of a(t) and a (t), a 
functional analysis is needed. However, the following method, "successive 
Fourier analysis," represents a practical way to accomplish the task. The 
first step is to Fourier-analyze the response over one period. For simplicity, 
a Fourier series with three terms will be used to explain the procedure of 
the modeling. Then 

C L = Aq + Aj cos0 + A 2 cos20 + A 3 cos30 

+ sin0 + B 2 sin20 + B 3 sin30 (4) 

The second step is to split the result into the form of eq. (3) with F(a,a ) 
and G(a,a ) being Fourier-analyzed again. The result after "successive 
Fourier Analysis" becomes 
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C L = A 0 + { CC[0,0] + CC[l,0]cc + CC[2,0]a 2 

+ DC[0,l]d + DC[l,l]aa + CC[0,2]d 2 }a 
+ { CS[0,0] + CS[l,0]a + CS[2,0]a 2 

+ DS[0,l]d + DS[l,l]aa + CS[0,2]d 2 }d (5) 

The detailed procedure of "successive Fourier analysis" is shown in 
Appendix 1. By collecting the same order terms of a, a and their products 
together, the result of becomes 

C L = A 0 + { CC[0,0]a + CS[0,0]d } 

+ { CC[l,0]a 2 + DC[0,l]act + CS[l,0]a<* 

+ DS[0,1] dc 2 } 

+ { CC[2,0]a 3 + DC[l,l]a 2 & 

+ CC[0,2]aa 2 + CS[2,0]a 2 d 

+ DS[l,l]acx 2 + CS[0,2]d 3 } (6) 

It can be seen that for each different frequency k with the same 
amplitude, there will be different response C L and different coefficients CC, 
DC, CS, and DS. To have practical applications, a general representation 
of these coefficients as a function of reduced frequency at a constant 
amplitude is needed. From the classical potential theory, it has been found 
that Pad6 approximants provide an accurate approximation of the 
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theoretical phase function. Therefore, Pad6 approximants will be used for 
the present model as phase functions to represent coefficients CC, DC, CS, 
and DS. Equation (6) is a useful form for determining stability derivatives 
based on forced oscillation tests. For applications to a maneuvering 
aircraft, the following representation of aerodynamic response based on the 
generalized indicial lift concept is more convenient. 

It is recalled that in the classical airfoil theory the circulatory lift is 
written as the product of Theodorsen’s circulation function and the quasi- 
steady lift. In the present nonlinear theory, the same form will also be 
adopted. For this purpose, eq. (4) (or the experimental oscillatory results) 
is rewritten in a complex form, as follows: 

C L = A 0 + (Aj - iBj) e ikt ’ + (A^* - iB 2 ) e i2kt ’ 

+ (A 3 - iB 3 ) e i3kt ’ (7) 

It should be kept in mind that only the real part of the response has 
physical meaning. The reason to put in the complex form is to benefit from 
the great mathematical convenience of the e ikt notation. If a is rewritten 


then eqs. (6) and (7) and the classical airfoil theory suggest that the 
response could be put in the following form involving the products of 
amplitude functions and phase functions as 


C L - C 0 (k) 

+ E X1 & + E 21 & + q * (H 31 a + h 21 6l ) * (1 - PD X ) 

+ E 12 & - E 22 & + C 2 * (H 12 a 2 + H 22 a & + H 32 d 2 ) 

* (i - pd 2 ) 

+ E 13 d + E 23 a + C 3 * (H^a 3 + H 23 a 2 ct + H 33 ad 2 + H i3 a 3 ) 

* (l - PD 3 ) 


(8) 


where PD is a Pad6 approximant with order 2 and is defined as 


P 3j (ik) 2 + P 2j (ik) 
P 3j (ik) 2 + (ik) + P tj 


E 11 a + E 21 & is the virtual-mass effect and accounts for the 
noncirculatory lift (ref. 2). In addition, H 21 , H 22 , H 23 , etc., are related to 
the pitch-rate effect. It should be noted that those terms inside the 
parentheses following C 1? C 2 , C 3 , such as (H^a + H 21 a ), represent the 
quasi-steady response and (1 - PDj) represents the unsteady aerodynamic 
lag in response. Therefore, the present assumed form for aerodynamic 
modeling encompasses the classical linear theory. 
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Cj are the reference values used to normalize the lift given by 
Aj - i Bj in the least squared-error method. A good choice for Cj is to use 
the same coefficient in the term as in the steady condition. Therefore, 
j is the index consistent with the exponent of the exponential term in eq. 
(7). For example if the j’s term in eq. (8) represents the coefficient of e^ , 
then j is 1. If the j’s term in eq. (8) represents the coefficient of e i2 ^ then 
j is 2, etc. The first term, C 0 (k), in eq. (8) is a constant term, supposedly 
a function of frequency. From available experimental data (ref. 19), it is 
found that an averaged constant can be used to represent C 0 (k) term as 
shown in Figure 1 for a delta wing. The unknown coefficients P^, P 2 j, P 3 j 
and P 4 j are calculated from the least squared-error method. E-q, E 21 , H^, 
H 12 , etc., will be obtained separately by minimizing the sum of squares of 
errors. This is equivalent to a two-level optimization method to determine 
the unknowns in eq.(8). That is, E, H, etc., are assumed first. Then P^, 
etc., are determined by minimizing the sum of squared errors. The values 
of E 11? H-q, etc., are varied next so that the sum of squared errors is 
minimized. It was found that this approach is more effective in 
determining a global minimum solution for the unknowns than a 
straightforward optimization (one level) method because of nonlinearity in 
the unknowns in the optimization problem. It should be noted that in the 
literature the phase function has been typically determined by the response 
to plunging motions. Therefore, those terms associated with a in eq. (8) do 
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not appear. This would very much simplify the mathematics of determining 
the Pad6 approximants. The details of the present method are discussed 
in the following. 


2.2 Least-Square Method 

By choosing proper values of E^, H^, H 12 , etc., in eq. (8), the 
corresponding Aj - i Bj term in eq. (7) is then divided by the amplitude 
function. The result will appear as 


v. 




( amplitude function) j 

Pjj (i-fc ) 2 + Pjj 
P 2j (ik ) 2 + (iie) + P 4j 


(9) 


If both sides of eq. (9) are multiplied by the denominator of the Pad6 
approximant and separated into real and imaginary parts, then 


Re s Pyk 2 - PyVjk 2 + P 4j V, - Wj k = 0 


(10a) 


and 


Im = P 2j k + P 3j Wjk 2 - P 4jWj - Vjk = 0 


(10b) 
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The sum of squared errors is defined as 


Err ■ £ Re{k ± ) 2 + £ ImikJ 2 ( 11 ) 


By equating the first derivatives of squared errors (eq. 11) with respect to 
variables P-jj, P 2 j, Pgj and P 4 j to zero, the unknown coefficients P^, P 2 j, Pgj 
and P 4 j can be determined by 


o E v ± k i 
o E*' E w ± k l 
-E v i k * E w ± k l E (vl k l+» 2 ikt) 
. E ^*1 -E E ( v * k l + > 


E v i k ± 

E 

E (V^kl + Wlkl) 

E 




E w i k i 

P 21 


E v i k * 



0 

P V. 


0 


( 12 ) 

where i varies over the range of input frequencies, and the mode subscript 
j on V and W has been omitted. 

2.3 Gradient Method 

After the unknown coefficients P^, P 2 j, P 3 j and P 4 j have been found, a 
one-dimensional gradient method is used to find E and H values which will 
make the sum of the squared errors minimum . The E or H value is 
perturbed first by a small amount AE or AH to find the gradient of the sum 
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of squared errors. If the gradient tends to reduce the error, then the E or 
H value is perturbed further until several iterations has been reached (it 
is set to be 5 iterations in the current program). After that, the same 
procedure is applied to other E or H. Then the whole procedure is repeated 
again until several iterations has been reached (it is set to be 10 in the 
current program). 

Finally, the response of C L is written as 


+ E X1 d + E 21 & + q * ( H lx a + H 2x i St) * (1 - PD X ) 

+ E 12 d + E 22 cl + C 2 * ( H 12 a 2 + H 22 ad + i/ 32 d 2 ) * (l - PD 2 ) 


+ q 3 d + e 23 & 


+ C 3 * ( H lz a 3 + H 23 a 2 d + H 32 ad 2 + H 43 i 3 ) * (l - PD 3 ) (13) 


It is easily seen that each term in the above equation is a product of an 
amplitude function and a phase function. The procedure to put oscillating 
response data into the form of eq. (13) is summarized in the next section. 

2.4 Summary of Numerical Procedure 

Step 1. Steady-state response analysis: 

Use Fourier analysis to analyze the steady-flow response over 
one period which is the same period as in the harmonic motion. Since in 
the steady flow the only variable in response is a, Fourier cosine series are 
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used. Then use eq. (A.3) to decompose the response into a polynomial of 
cos(kt’) or a. If the final result for the steady response is written as 

C L = I 0 + I^a + I 2 a 2 + (14) 

then the coefficients C 0 , and Cj in eq. (14) are the same value as I Q , and Ij 
where j = 1,2,3. 

Step 2. Unsteady-response analysis: 

Use Fourier analysis to analyze the harmonic motion response 
for different frequencies over one period. For each frequency, the response 
should be in the same form as in eq. (7). 


Step 3. Constant-term analysis: 

From step 2, if constant terms appear in the Fourier analysis 
then C ave is calculated from each constant term due to different frequencies 
as 


E 

n*l 


( 15 ) 


16 



where N is the number of frequencies used in step 2 and A^(k n ) are the 
constant terms due to different frequencies in Fourier analysis. 


Step 4. Amplitude-phase identification: 

In this step, the coefficients Aj and Bj calculated from step 2 are put 
into the form of a product of amplitude functions and phase functions as in 
eq.(14). The procedure is as follows: 

4-a Set the initial guess for E or H values. 

4-b Use the least-square method to find unknown coefficients 

P lj> P 2j * P 3j and p 4j- 

4-c Use gradient method to find better E or H values. 

Repeat steps 4-a to 4-c until the sum of square errors reached 
minimum or the setting iteration limit has been reached. 

Although three-term Fourier series are used in the above, the 
procedure is applicable to any number of Fourier terms. 


2.5 Indicial Formulation 

In linear theory, the reciprocal relations (or Fourier summation) has 
been used to calculate the indicial response. However, in nonlinear theory 
those relations can not be applied. As Tobak (ref. 7) mentioned in his 
paper, the aerodynamic response due to a step change should reach steady- 
state value asymptotically at subsonic speeds. In linear theory, these 
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asymptotic relations are represented by exponential functions (refs. 2 and 
5), and these exponential functions are calculated through the phase 
function. Therefore, the phase function in the current nonlinear modeling 
will be converted into exponential function in time domain but keep the 
amplitude function unchanged. If eq. (13) is rewritten for m-terms Fourier 
series as 

c - c 

+ E lx d + E 21 &+ C x * (H X1 a + H 21 6l) * (1 - PD X ) 

+ E 12 & + E 22 & + C 2 * (tf 12 a 2 + H 22 Cl6l + H 22 d 2 ) * (1 -PD 2 ) 

+ E 12 6l + E 23 & + C 3 * (tf 13 a 3 + H 22 cl 2 6l + H 32 ai 2 + tf 43 d 3 ) 

* (1-PD 3 ) (16) 

+ 

m 

+ E + E 2^ 

>1 

+ ( amplitude function) j * ( phase function) j 


then by Fourier inversion of the Pad6 approximants, the integrand in eq. 
(2) can be obtained from the following expression: 
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L ‘indlciMl 


q *(H 11 a + H 21 ft) * (1 - a 31 e'“ Jl 
+ c 2 * (H 12 ff 2 + tf 22 ad + q 2 d 2 ) 


3 21 S’' 41 ' ) 


-a 13 2 c 


* (1 - a 12 e' a,J 


- a 22 e 




) 


y; (amplitude function) j 
£ i 


* (l - a XJ e' a ^ Jt/ - 


a 2J e~ a * jjt ') 


( 17 ) 


where the coefficients a^, a 2 j, a 3 j and a 4 j, are calculated from Pade 
approximants ascorresponding indicial response for nonlinear theory is 
defined as 


Pjj (ik ) 2 + Pzj 
P 3j (ik ) 2 + ik + P Aj 


** a U + a 2 J 

iJc + a 3j iic + a 4 _^ 

i ( jk) a 3j i (j'Jc) a 2i 

i ( j£) + ja 3i i (j'Jc) + a 4 ^ 


Then the generalized response for arbitrary motion is obtained by time 
integration of eq. (2). 
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Chapter 3 

Results and Discussion 


3.1 Linear Results 

Several cases in the two-dimensional and three-dimensional linear flow 
have been studied to verify the proposed method of aerodynamic modeling. 

3.1.1 Two-Dimensional Flow 

The first case studied is a 2-D flat plate oscillating in the incompressible 
flow. The amplitude is 57.3 degree (one radian) in angle of attack for the 
airfoil. Therefore it oscillates from 57.3 degree of angle of attack to -57.3 
degree of angle of attack then back to 57.3 degree of angle of attack for one 
cycle with respect to midchord, i.e. 

dj = 0.0 + 1.0 cos(kt’) (in radian) 

The steady lift is already known from the linear theory (ref. 2) as 2m, 
and the oscillating complex lift is taken from a 2-D unsteady QVLM 
program (ref. 19) as input data for the current model and for comparison. 
Through numerical experimentation, it is found that only six frequencies 
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are needed to have accurate results. Through the modeling procedure as 
summarized in section 2.4, the lift can be written as 


c { = E n a + E 21 & + ( H-qCc + H 21 a ) (1.0 - PD 1 ) 


where 


PD, 


P^Uk ) 2 P 21 (iJc) 
P 31 (iJc ) 2 + ik + P 41 


and the values of coefficients E n , E 21 , H 11? H 21 and P i;L are listed in Table 
1. The results for the lift coefficients are plotted in Figure 2 for different 
numbers of frequencies used. Compared with the aerodynamic responses 
by the 2-D UQVLM program, the numerical results from modeling show 
excellent agreement. 

Two Mach numbers, 0.2 and 0.4, are chosen in the 2-D compressible 
flow to verify the current model. A two-dimensional unsteady QVLM (ref. 
19) program is again used to calculate the complex lift as input data for the 
current model and for comparison. The same frequencies as in the 
incompressible flow are used as input also. The results for coefficients E-q, 
E 21 , H^, H 21 and P^ are listed in Table 1. The aerodynamic responses c^ 
calculated by the model are plotted in Figures 3 and 4 to compare the 
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results from 2-D unsteady QVLM. These figures show that the numerical 
results by modeling are very accurate. 

3.1.2 Three-Dimensional Flow 

The same Mach numbers (includes 0) in the 2-D flow are used in 3-D 
attached flow to verify the current model. The geometry is a 70 degree 
delta wing which oscillate from zero degree angle of attack to twenty degree 
angle of attack i.e. 

= 0.1745329 + 0.1745329 cos kt’ (in radian) 

This means that the mean angle of attack is ten degree (0.1745329 radian) 
and the amplitude of the oscillation is ten degree (0.1745329 radian). The 
input aerodynamic responses are calculated from a 3-D unsteady QVLM 
program (ref. 20). In the program, the total lift is the sum of steady lift at 
the mean angle-of-attack plus unsteady lift. Since the steady lift is the 
same for every term, only the unsteady lift is used in the modeling and for 
comparison. Through numerical experimentation, it is found that the 
responses at low frequencies do not change significantly, which results in 
inaccurate modeling. To have accurate approximation, high frequencies’ 
responses are needed. Seven reduced frequencies (k = 0.01, 0.1, 0.2, 0.6, 
1.0, 2.0, 2.5) are used as input data in the 3-D attached flow cases. The 
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results for the coefficients of the modeling are listed in Table 1. The 
responses C L from modeling are plotted in Figures 5 to 7 to compare with 
the results from 3-D unsteady QVLM program. All of these figures show 
very good agreement. 

3.2 Nonlinear Result 

The experimental data (ref. 21) for a 70-degree delta wing in pitching 
oscillation is used to validate the current aerodynamic model. The angle 
of attack which describes the pitching motion is given as 
= 27.5 - 27.5 cos kt’ (in degree) 

which means the delta wing oscillates from zero degree angle of attack to 
55-degree angle of attack then back to zero degree angle of attack for one 
cycle. The reduced frequency k is nondimensionalized based on wing’s root 
chord. Five sets of data corresponding to five different frequencies are 
available and they will be used as the input data to calculate the 
coefficients for the current aerodynamic model. Five terms in the Fourier 
series are used for the current aerodynamic model. The calculated 
coefficients for the current aerodynamic model are listed in Table 2, and the 
response Cl is written as 
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C avg 

+ E 12 d + e 21 & + q* (H n a + h 21 &) * (l - PD X ) 

+ £12* + *22* + ^12« 2 + #22®* + ^32* 2 ) * <1 “ ^ 

+ E 12 6 : + E 22 & + C 3 *(H 13 a 3 + H 22 a 2 d + tf 33 ad 2 + tf 43 d 3 ) 

* (1 - PD 3 ) 

+ i? 14 d + E 2i & + C 4 * (tf 14 a 4 + tf 24 a 3 d + f/ 34 a 2 d 2 + tf 44 ad 3 + tf 54 d 4 

* (1 - pd 4 ) 

+ q 5 d + E 25 & 

+ C 5 * {h 15 a 5 + tf 2 5 a 4 a + H 35 a 3 a 2 + f/ 45 a 2 d 3 + ff 55 ad 4 + tf 65 d 5 ) 

* (1 - P.D 5 ) 

where PDj are the Pade approximants. 

The result from modeling is plotted in Figure 8, which shows reasonable 
agreement with experimental data for each frequency. 

3.2.1 Indicial Formulation 

Note that eq. (2) is valid for arbitrary motion. To check its validity in 
the nonlinear theory, two oscillatory cases in the last section will be used. 
That is, by assuming oscillatory motion in eq. (2), the time-integrated lift 
response should agree with the forced-oscillation data. As indicated 
earlier,, the integrand of eq. (2) can be written as follows: 
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= C 1 * (if H a + H Z1 &) * (1 - eX Pl> 

+ C 2 * {H 12 a 2 + H 22 ad + ff 32 d 2 ) * (1 - exp 2 ) 

+ C 3 * (if 13 a 3 + H 23 a 2 d + ff 33 ad 2 + tf 43 d 3 ) * (l-exp 3 ) 

+ C 4 * (/f 14 a 4 + H 2i a 3 d + H 3i a 2 a 2 + i/ 44 ad 3 + tf 54 d 4 ) 

* (1 - exp 4 ) 

+ C 5 * (£T 15 a 5 + H 25 a*a + H 35 a 3 a 2 + /f 45 a 2 a 3 + ff 55 od 4 + tf 65 d 5 ) 

* (l - exp 5 ) 


where the functions expj are converted from Pad6 approximants as 


exp j — &ij 




* 2 ] 




By differentiating with a and a , C T , Ct are obtained as follows: 

a a 

c l % ~ c i * #11 * d ~ exp x ) 

+ C 2 * (2 tf 12 a + ) * (1 - exp 2 ) 

+ C 3 * (3H 13 a 2 + 2H 23 od + H 33 <t 2 ) 

* (1 - exp 3 ) 

+ C 4 * (4if 14 a 3 + 3ff 24 a 2 d + 2 i/ 34 ad 2 + ff 44 d 3 ) 

* (1 - exp 4 ) 

+ C s * (5tf 15 a 4 + 4tf 25 a 3 d + 3tf 35 a 2 d 2 + 2/f 45 ad 3 + tf 55 d 4 ) 

* (1 - exp 5 ) 
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= q * ff 21 * (1 - exPi) 

+ c 2 * (W 2 2 a + 2tf 32 d) * (1 - exp 2 ) 

+ C 3 * (if 23 a 2 + 2H J3 aa + 3if 43 d 2 ) * (1 - exp 3 ) 

+ C 4 * (tf 24 a 3 + 2 tf 34 a 2 d + 3 tf 44 ad 2 + 4tf 54 d 3 ) 

* (1 - exp 4 ) 

+ c 5 * (tf 25 a 4 + 2tf 35 a 3 a + 3tf 45 a 2 d 2 + 4tf 55 ad 3 + 5tf 65 d 5 ) 

* (l - exp 5 ) 


The function C L (0) in eq. (2) includes the initial conditions and 
virtual mass effect. In the present case, C L (0) is calculated as 

C L (0) = C L (t,a(0),d (0)) + C ave 

+ E^a + E 21 & + E 12 a 2 + E 22 d2 + E 13 a 3 + E 23 & 3 
+ E 14 a 4 + E 24 ix 4 + E 15 a 5 + E 25 & 5 

where the subscript to the a-terms indicates the order of the harmonics. 
For example, c ^ is proportional to the second harmonics. Simpson’s 1/3- 
rule is used to integrate eq. (2). Since the angle of attack is set to be a 
complex number (cos(kt’)+i sin(kt’)) in oscillating cases, only the real part 
of the integrated lift is taken. The lift by integrating eq. (2) for a 70-deg. 
oscillating delta wing with frequencies k=0.098 and k=0.165 are plotted in 
Figure 9. Compared with the lift from aerodynamic modeling, the 
integrated lift shows good agreement. 
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RESULTS FROM FOURIER ANALYSIS 
NUMERICAL APPROXIMATION 




LIFT COEFFICIENT 



REDUCED FREQUENCY k 

FIG. 2 COMPARISON OF NUMERICALLY APPROXIMATED COMPLEX 
LIFT COEFFICIENT WITH 2-D UNSTEADY QVLM 
RESULTS AT M-0.0 
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LIFT COEFFICIENT 



FIG. 3 COMPARISON OF NUMERICALLY APPROXIMATED COMPLEX 
UFT COEFFICIENT WITH 2-0 UNSTEADY QVLM 
RESULTS AT M-0.2 
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LIFT COEFFICIENT 



REDUCED FREQUENCY k 

FIG. 4 COMPARISON OF NUMERICALLY APPROXIMATED COMPLEX 
UFT COEFFICIENT WITH 2-D UNSTEADY QVLM 
RESULTS AT M-0.4 
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LIFT COEFFICIENT 
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M - 0.0 (7 INPUT FREQUENCIES) 

0 RESULTS FROM 3-D UQVLM (in-phase) 

A RESULTS FROM 3-0 UQVLM (out-of-phase) 

NUMERICAL APPROXIMATION 
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FIG. 5 COMPARISON OF NUMERICALLY APPROXIMATED COMPLEX 
LIFT COEFFICIENT WITH 3-0 UNSTEADY QVLM 
RESULTS AT M-0.0 
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LIFT COEFFICIENT 



0.0 0.3 1.0 1.3 2.0 2.3 

REDUCED FREQUENCY k 

FIG. 6 COMPARISON OF NUMERICALLY APPROXIMATED COMPLEX 
LIFT COEFFICIENT WITH 3-D UNSTEADY QVLM 
RESULTS AT M-0.2 
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LIFT COEFFICIENT 



REDUCED FREQUENCY k 

FIG. 7 COMPARISON OF NUMERICALLY APPROXIMATED COMPLEX 
UFT COEFFICIENT WITH 3-D UNSTEADY QVLM 
RESULTS AT M-0.4 
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LIFT COEFFICIENT LIFT COEFFICIENT 




ANGLE OF ATTACK 

FIG. 6 COMPARISON OF NUMERICALLY APPROXIMATED LIFT 
COEFFICIENT WITH EXPERIMENTAL DATA FOR 
AN OSCILLATING 70-DEG. DELTA WING 
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ANGLE OF ATTACK 
RG. 8 CONTINUED 
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UFT COEFFICIENT LIFT COEFFICIENT 




ANCLE OF ATTACK 
FIG. 8 CONCLUSION 
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LIFT COEFFICIENT LIFT COEFFICIENT 




NONDIMEN SIONAL TIME ? 

FIG. 9 COMPARISON OF LIFT COEFFICIENT CALCULATED FROM 
NUMERICAL APPROXIMATED MODELING AND FROM 
INDICIAL LIFT TIME INTEGRATION FOR AN 
OSCILLATING 70-DEG DELTA WING 
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Appendix 

Successive Fourier Analysis 


The first step of "successive Fourier analysis" is to Fourier- analyze the 
response over one period. For simplicity, a Fourier series with three terms 
will be used to explain the procedure of the modeling. Then 


C L = A 0 + A i cosG + A 2 cos29 + A 3 cos30 

+ B^^ sin0 + B 2 sin20 + B 3 sin30 (Al) 

where 

a = "^r/o Cl ** 

A„ = — [ 2n C L cos (nQ) d0 

n n Jo L 

B. - — f 2n C L sin (n0) dB 

n 71 Jo L 

n=l,2,3 and 0=Jct / 


( A2 ) 


Once the coefficients of Aq, A^ and B n have been found, the next step is 
to split the coefficients into two groups by using the following formulas, 


A.l 



cos nd = C(n,0) cos "6 - C(n, 2) cos n_2 0 sin 2 0 
+ C(n, 4) cos n_4 0 sin 4 0 + 

sin n0 = C(n,l) cos n_1 0 sin0 - C{n, 3) cos n_3 0 sin 3 0 

- C(n, 5) cos"' 5 0 sin 5 0 + (A2) 


where 


C{n, m) 


n! 

( n-m ) ! m\ 


and 


n\ =1*2*3*4** ****n 


Therefore, the response of Cl becomes 

C L = A 0 + A 1 cos0 + A 2 [cos 2 0 - sin 2 0] 

+ A 3 [cos 3 0 - 3cos0 sin 2 0] 

+ B-^ sin0 + B 2 [2 cos 0 sin0] 

+ B 3 [3cos 2 0 sin0 - sin 3 0] 

= Aq + [A-l + A 2 cos0 + A 3 (cos 2 0 - 3sin 2 0)] cos0 
+ [B1 + 2B 2 cos0 + B 3 [3cos 2 0 - sin 2 0] sin0 
= Aq + F(cos0 , sin0) cos0 + G(cos0 , sin0) sin0 

Perform the Fourier analysis again for functions F(cos0 , sin0) and 
G(cos0 , sin0) by using Fornier series with the same terms as in the first 
step. Then 


A.2 



F(cos0,sin0) = Fq + FA 1 cos0 + FA 2 cos20 + FA 3 cos30 
+ FBj^ sin0 + FB 2 sin20 + FB 3 sin30 
G(cos0,sin0) = GO + GA 1 cos0 + GA 2 cos20 + GA 3 cos30 
+ GBj^ sin0 + GB 2 sin20 + GB g sin30 

where 


F n = 

i 

r f d e 

Gn = 

1 

f 2 * G d< 

u 

271 

Jo 


0 

271 

Jo 

3 

ii 

— / 
n J c 

2 n 

1 

F d 0 

GA n = 

vJ 

r Gde 
0 

FB n = 

— / 
71 J C 

271 

1 

F d 8 

GB n = 

iJ 

r Gde 
0 

n = 

1,2, 

3 






Using eq. (A3) again, then 

F(cos0,sin0) = F Q + FA 1 cos0 + FA 2 [cos 2 0 - sin 2 0] 
+ FA 3 [cos 3 0 - 3 cos0 sin 2 0] 

+ FB 1 sin0 + FB 2 [2 cos 0 sin0] 

+ FB 3 [3cos 2 0 sin0 - sin 3 0] 

G(cos0,sin0) = G 0 + GAj cos0 + GA 2 [cos 2 0 - sin 2 0] 
+ GA 3 [cos 3 0 - 3 cos0 sin 2 0] 

+ GBj^ sin0 + GB 2 [2 cos 0 sin0] 

+ GB 3 [3cos 2 0 sin0 - sin 3 0] 


A.3 


Therefore 


Cl = A 0 + {Fq + FA-^ cosO + FA 2 [cos 2 0 - sin 2 0] 

+ FA 3 [cos 3 0 - 3 cos0 sin 2 0] 

+ FB 2 sin0 + FB 2 [2 cos 0 sin0] 

+ FB 3 [3cos 2 0 sin0 - sin 3 0]} cos0 
+ {Gq + GA-^ cos0 + GA 2 [cos 2 0 - sin 2 0] 

+ GA 3 [cos 3 0 - 3 cos0 sin 2 0] 

+ GB]^ sin0 + GB 2 [2cos0 sin6] 

+ GB 3 [3cos 2 0 sin0 - sin 3 0]} sin0 

All the terms associated with cos n 0 on the right hand side of the above 
equation are divided by (a 0 ) n and the terms associated with sin n 0 are 
divided by (-ka 0 ) n . After rearrangement, the response of Cl becomes 

C L = A 0 + {CC[0,0] + CC[1,0] a + CC[2,0] a 2 + CC[3,0]a 3 
+ DC[0,1] d + DC[1,1] ad + DC[2,1] a 2 d 
+ CC[0,2] d 2 + CC[1,2] ad 2 + DC[0,3] d 3 } a 
+ {CS[0,0] + CS[1,0] a + CS[2,0] a 2 + CS[3,0]a 3 
+ DS[0,1] d + DS[1,1] ad + DS[2,1] a 2 d 
+ CS[0,2] d 2 + CS[1,2] ad 2 + DS[0,3] d 3 } a (A4) 


A.4 


where 


CC[n,m ] 
CS[n,/n] 


[«0 n (-iea®) ] 
[«o a ( -kan) ] 


DC[n,m ] 
DS[n,m] 


ggn*m 

[o 0 n (-Jca 0 ")l 




[a<? (-ica®) ] 


and the coefficients CC[n,m], DC[n,m], CS[n,m] and DS[n,m] are zeros for 
n+m t 3. Comparing with eq. (Al), it is obtained that 

F 0 = A 0 

F(a, d) = CC[0,0] + CC[1,0] a + CC[2,0] a 2 

+ DC[0,1] cx + DC[1,1] act + CC[0,2] a 2 
G(a, fit) = CS[0,0] + CS[1,0] a + CS[2,0] a 2 

+ DS[0,1] a + DS[1,1] act + CS[0,2] d 2 

Finally, collecting the same order terms together, then 
C L = Ao 

+{ CC[0,0]a + CS[0,0]ct } 

+{ CC[l,0]a 2 + DC[0,l]ad +CS[l,0]ad + DS[0,1] d 2 } 

+{ CC[2,0]a 3 +DC[l,l]a 2 d + CC[0,2]otd 2 +CS[2,0]a 2 d 
+DS[l,l]ad 2 + CS[0,2]d 3 } (A5) 


A.5 



